%
%
%

clear all

dynare olg_ann.mod nostrict

%%

n_ = M_.endo_nbr ; % Number of endogenous variables
l_ = M_.exo_nbr ;  % Number of exogenous variables
n_ = n_+1 ;        % Constant

SET.variable.n_ = n_ ; 
SET.variable.l_ = l_ ; 
SET.nparam      = M_.param_nbr ;

var_names   = cellstr(M_.endo_names) ;
param_names = cellstr(M_.param_names) ;
exo_var_names = cellstr(M_.exo_names) ;

for ii=1:n_-1
    eval(['SET.variable.',var_names{ii},'=',int2str(ii),';']) ;
end

for ii=1:l_
    eval(['SET.variable.shock.',exo_var_names{ii},'=',int2str(ii),';']) ;
end

for ii=1:M_.param_nbr
    eval(['SET.param_names.',param_names{ii},'=',int2str(ii),';']) ;
end

ky = oo_.steady_state(SET.variable.k) ./ oo_.steady_state(SET.variable.y) ;
cy = oo_.steady_state(SET.variable.c) ./ oo_.steady_state(SET.variable.y) ;
ln = oo_.steady_state(SET.variable.l) ;
wly = oo_.steady_state(SET.variable.w) * oo_.steady_state(SET.variable.l) * oo_.steady_state(SET.variable.n) / oo_.steady_state(SET.variable.y) ;
r = 100*(-1+oo_.steady_state(SET.variable.R)) ;

fprintf(sprintf(...
    ' KoY  .... %1.3f \n CoY  .... %1.3f \n LoN  .... %1.3f \n wLoY .... %1.3f \n R    .... %1.3f \n', ...
    ky, cy, ln, wly, r)) ;

%% Lifecycle figures

for tt=1:203
for jj=1:N_g
    eval(['c_profile(',int2str(jj),',' int2str(tt), ') = c', int2str(jj),'(',int2str(tt),') ;']) ;
end

for jj=1:N_g-1
    eval(['k_profile(',int2str(jj),',' int2str(tt), ') = k', int2str(jj),'(',int2str(tt),') ;']) ;
end
end

for tt=1:203
for jj=1:N_g
    eval(['n_profile(',int2str(jj),',' int2str(tt), ') = n', int2str(jj),'(',int2str(tt),')  ;']) ;
    eval(['n_profile_(',int2str(jj),',' int2str(tt), ') = n', int2str(jj),'_(',int2str(tt),') ;']) ;
end
end

for tt=1:203
for jj=1:N_g-1
    eval(['wealth_dist(',int2str(jj),',' int2str(tt), ') = n', int2str(jj),'_(',int2str(tt),')*k', int2str(jj),'(',int2str(tt),')  ;']) ;
    eval(['income_dist(',int2str(jj),',' int2str(tt), ') = n', int2str(jj),'(',int2str(tt),')*n', int2str(jj),'_(',int2str(tt),')*w(',int2str(tt),')*a', int2str(jj), ' ;']) ;
end
end

%% Aggregate shocks

if length(vvec_end)<203
    for ii=132:203
        vvec_end = [vvec_end vvec_end(:,131)] ;
        gam_mat = [gam_mat gam_mat(:,131)] ;
    end
end

for tt=1:203
    for jj=1:80
        phi_js(jj,tt) = 1;
    end
end

for tt=1:203
    B_s(tt) = sum(n_profile_(:,tt).*(avec'.*(n_profile(:,1)>0)))  ;
end

for tt=1:203
    avec_norm(:,tt) = avec ./ B_s(tt) ;
end

for tt=1:203
    A_s(tt) = sum(n_profile_(:,tt).*(avec_norm(:,tt).*(n_profile(:,1)>0)).^(1+1/varphi).*(phi_js(:,tt).*vvec_end(:,1)).^(-1/varphi)) ;
    A_s(tt) = A_s(tt) / sum(n_profile_(:,tt).*(avec_norm(:,tt).*(n_profile(:,1)>0)).^(1/varphi).*(phi_js(:,tt).*vvec_end(:,1)).^(-1/varphi)) ;
end

for tt=1:203
    phi_s(tt) = (sum(n_profile_(:,tt) .* phi_js(:,tt).^(1/nu)))^nu ; 
end

for tt=1:203
    v_s(tt) = (sum(n_profile_(:,tt).*(avec_norm(:,tt).*(n_profile(:,1)>0)).^(1+1/varphi).*(phi_js(:,tt).*vvec_end(:,1)).^(-1/varphi)))^(-varphi) ;
end

tau = oo_.endo_simul(SET.variable.tau,:) ;

for tt=1:203
    lgv_vec(tt) = mean(gam_mat(:,tt).*gam_mat(:,tt)) ;
end

timevec_le = 1860:2070 ;
life_exp2 = life_exp(find(timevec_le==1940):find(timevec_le==2070)) ;
life_exp2_pr = 1 ./ (life_exp2*1) ;
life_exp2_pr_ch = ([0 (life_exp2_pr(2:end)) - (life_exp2_pr(1))]) ;
life_exp2_pr_ch = life_exp2_pr_ch+1 ;

lgv_vec = life_exp2_pr_ch ;

%tauv = tau ;

% save demographic adjustment factors for matlab here
%save adj_factors A_s v_s phi_s B_s tauv lgv_vec
   
%% Data 

% From MFP, in data/'BLS capital output ratio' (Private Non-Farm Business
% Sector)
capout=[2.136
2.213
2.089
2.047
2.058
2.026
2.109
2.016
2.063
2.090
2.186
2.065
2.093
2.106
2.035
2.014
1.962
1.924
1.904
1.984
1.971
2.012
2.112
2.119
2.073
2.040
2.183
2.307
2.223
2.186
2.146
2.194
2.340
2.415
2.607
2.538
2.460
2.486
2.507
2.519
2.500
2.507
2.554
2.644
2.603
2.603
2.576
2.597
2.603
2.609
2.630
2.653
2.704
2.808
2.848
2.833
2.778
2.758
2.763
2.781
2.890
3.038
2.955
2.919
2.869
2.870
2.842
2.836];

r_data = (1/3) * (capout.^(-1)) ;

capout_model = oo_.endo_simul(SET.variable.k,:)./oo_.endo_simul(SET.variable.y,:) ;
caplab_model = oo_.endo_simul(SET.variable.k,:)./oo_.endo_simul(SET.variable.n,:) ;

emp_pop_ratio_dates = 1948:.25:2020.75;

% From FRED 
% Units:  Percent, Seasonally Adjusted
% Frequency:  Monthly
% The series comes from the 'Current Population Survey (Household Survey)'
% The source code is: LNS12300000
% Suggested Citation:
% U.S. Bureau of Labor Statistics, Employment-Population Ratio [EMRATIO]
emp_pop_ratio = [56.5
56.6
56.8
56.6
56.1
55.4
55.1
55.3
55.1
55.9
56.5
56.8
57.2
57.3
57.4
57.4
57.5
57.2
57.1
57.3
58.0
57.3
57.1
56.3
55.9
55.4
55.2
55.4
55.7
56.3
57.1
57.4
57.5
57.5
57.6
57.4
57.4
57.2
57.1
56.6
55.6
55.3
55.3
55.5
55.7
56.3
56.1
56.0
55.9
56.4
56.2
55.9
55.6
55.3
55.2
55.4
55.6
55.6
55.6
55.3
55.2
55.4
55.4
55.4
55.5
55.9
55.7
55.6
55.8
56.1
56.3
56.5
56.6
56.8
57.0
57.3
57.0
57.1
57.4
57.5
57.2
57.7
57.5
57.6
57.8
57.9
58.1
58.1
57.9
57.6
57.2
56.9
56.6
56.5
56.6
56.7
56.8
57.0
57.0
57.2
57.5
57.8
57.9
58.2
58.2
58.0
57.8
57.3
56.2
55.9
56.1
56.1
56.5
56.9
57.0
57.0
57.2
57.8
58.0
58.5
58.8
59.3
59.4
59.8
60.0
59.8
59.9
60.0
59.9
59.1
58.8
59.0
59.2
59.4
59.0
58.5
58.2
58.0
57.7
57.3
57.1
57.5
58.2
58.6
59.0
59.6
59.7
59.8
60.0
60.0
60.1
60.4
60.5
60.6
60.8
60.9
61.1
61.4
61.7
61.9
62.0
62.2
62.4
62.6
62.9
62.9
63.0
63.0
63.2
63.0
62.7
62.3
61.9
61.8
61.6
61.4
61.4
61.5
61.5
61.4
61.4
61.7
61.8
61.9
62.2
62.4
62.5
63.0
63.1
62.8
62.8
62.8
62.9
63.1
63.3
63.4
63.5
63.7
63.9
64.0
64.0
64.1
64.0
64.2
64.3
64.2
64.2
64.4
64.6
64.5
64.2
64.3
64.3
63.8
63.5
63.0
62.8
62.8
62.8
62.5
62.5
62.3
62.1
62.2
62.3
62.3
62.4
62.4
62.4
62.7
62.8
62.8
63.0
63.1
63.1
63.3
63.3
63.0
62.8
62.8
62.8
62.5
62.0
61.4
60.3
59.6
59.0
58.5
58.5
58.6
58.5
58.3
58.4
58.3
58.3
58.5
58.5
58.5
58.5
58.7
58.6
58.6
58.7
58.5
58.8
58.9
59.0
59.3
59.2
59.4
59.3
59.4
59.8
59.7
59.8
59.7
60.0
60.1
60.2
60.1
60.3
60.5
60.4
60.5
60.7
60.6
60.8
61.0
60.7
52.9
56.1
57.4];

% Net Saving
% Suggested Citation:
% U.S. Bureau of Economic Analysis, Net saving as a percentage of gross 
% national income [W207RC1Q156SBEA], retrieved from FRED, Federal Reserve 
% Bank of St. Louis; https://fred.stlouisfed.org/series/W207RC1Q156SBEA, December 31, 2020.

% population growth data: SPPOPGROWUSA

% population growth projections: https://www.census.gov/data/tables/2017/demo/popproj/2017-summary-tables.html
% census

netsavingmodel = (k(2:end)-k(1:end-1))./y(1:end-1) ;

load ../data/demographic/netsavingdata
load ../data/demographic/popgrowthdata
load ../data/demographic/popgrowthprojdata

date_netsavingmodel=1941:2142;
netsavingmodel(find(date_netsavingmodel==1961):find(date_netsavingmodel==2019)) = ...
    netsavingmodel(find(date_netsavingmodel==1961):find(date_netsavingmodel==2019))+popgrowth;
netsavingmodel(find(date_netsavingmodel==1960))=netsavingmodel(find(date_netsavingmodel==1960)) + popgrowth(1);
netsavingmodel(find(date_netsavingmodel==2020):find(date_netsavingmodel==2060)) = ...
    netsavingmodel(find(date_netsavingmodel==2020):find(date_netsavingmodel==2060)) + popgrowthproj/100 ;

%% Figure 4: Macro Trends Due to Demographics

figure ;

set(gcf, 'DefaultAxesLineWidth', 1);
set(gcf, 'DefaultLineLineWidth', 2);
set(gcf, 'DefaultAxesTickLabelInterpreter','latex'); 
set(gcf, 'DefaultLegendInterpreter','latex');
set(gcf, 'DefaultAxesFontSize',12);
set(gcf, 'Position', [1.6353e+03 712 560*(2) 420*(2/3)])

subplot(1,3,1) ; hold on ;
plot(emp_pop_ratio_dates,emp_pop_ratio)
plot(1940:2142,l*100,'-.') 
xlim([1960 2060]) ;
title('A. Employment to Population, \%','Interpreter','latex')
h=legend('Data','Model','location','northeast');
set(h,'Interpreter','latex','FontSize',11) ;
ax = gca ;
ax.YGrid='on';
box('off') ;
xlim([1980 2040]) ;
ylim([50 70]) ;

subplot(1,3,2) ; hold on ;
plot(1948:2015,capout);
plot(1940:2142,capout_model,'-.') ;
xlim([1980 2040]) ;
title('B. Capital to Output','Interpreter','latex')
ax = gca ;
ax.YGrid='on';
box('off') ;

subplot(1,3,3) ; hold on ;
plot(date_netsaving,netsaving/100)
plot(1941:2142,netsavingmodel+0.0134,'-.') % Add estimated TFP growth
line([1940 2142],[0 0],'Color','k','LineWidth',1) ;
xlim([1980 2040]) ;
title('C. Net Saving to Output','Interpreter','latex')
ax = gca ;
ax.YGrid='on';
box('off') ;
